A GeoestatÃstica é uma ferramenta na análise espacial de dados, desempenhando um importante papel na agricultura moderna. Ela combina princÃpios estatÃsticos com informações geográficas para fornecer uma compreensão mais profunda da variabilidade espacial de várias variáveis ligadas ao solo, à s plantas e mesmo à atmosfera.
É definida como um ramo da estatÃstica aplicada que lida com a análise e interpretação de dados espaciais. Ela permite a modelagem da variabilidade espacial de fenômenos, descrevendo padrões espaciais, correlações e distribuições geográficas. No contexto agronômico, é uma ferramenta para compreender a heterogeneidade do solo e outros fatores que influenciam nos mais diversos fenômenos associados.
O objetivo deste curso é capacitar os alunos a utilizar os princÃpios da GeoestatÃstica para analisar e interpretar dados espaciais. Ao final do curso, os participantes serão capazes de aplicar métodos geoestatÃsticos para mapear variáveis do solo, identificar padrões de distribuição espacial e tomar decisões informadas sobre práticas agrÃcolas.
Considere os dois conjuntos de dados \(A\) e \(B\) e suas respectivas estatÃsticas básicas.
| EstatÃstica | A | B |
|---|---|---|
| n | \(15251\) | \(15251\) |
| média | \(100,0\) | \(100,0\) |
| desvio padrão | \(20,0\) | \(20,0\) |
| mediana | \(100,35\) | \(100,92\) |
| Percentil 10 | \(73,89\) | \(73,95\) |
| Percentil 90 | \(125,61\) | \(124,72\) |
De acordo com estas evidências os dois conjuntos de dados são bastante semelhantes.
Comparação de seus respectivos gráficos
de contornos (padrões espaciais).
O conjunto \(A\) é mais acidentado que o conjunto de dados \(B\).
Não se pode afirmar que o conjunto de dados \(A\) é mais variável do que o conjunto \(B\), haja visto que os desvios padrões dos dois conjuntos de dados foram iguais.
O conjunto \(A\) muda mais rapidamente no espaço do que o conjunto \(B\).
As zonas contÃnuas altos valores (em vermelho) e as zonas de baixos valores (em azul) são, em média, menores para o conjunto \(A\) do que do conjunto \(B\).
Estas diferenças podem ter um impacto significativo no delineamento amostral, na caracterização do local e na predição espacial em geral.
Portanto, não é surpresa que as estatÃsticas comuns falhem em reconhecer as diferenças da variável respostas dos dois exemplos de conjuntos de dados, isso porque a estatÃstica descritiva e os histogramas não incorporam a localização espacial dos dados em seus cálculos.
Podemos entender a geoestatÃstica como estudos de fenômenos que variam no espaço e/ou no tempo, ou seja, é uma coleção de técnicas numéricas, que lidam com a caracterização de atributos espaciais, permitindo a descrição da continuidade espacial de fenômenos naturais e fornece adaptações das técnicas da regressão para o entendimento desta continuidade (ISAAKS; SRIVASTAVA, 1989; GOOVAERTS, 1997).
Embora ela tenha sua origem na mineração, a geoestatÃstica é uma parte básica de muitas disciplinas cientÃficas incluindo geosfÃsica, ciências agronômicas, hidrologia e engenharia ambiental. A parte central da geoestatÃstica é a ideia de que medidas mais próximas tendem a serem mais parecidas do que valores observados em locais distantes. A geoestatÃstica fornece métodos para quantificar esta correlação espacial e incorporá-la na estimação e na inferência (GOTWAY; HARTFORD).
O variograma (ou semivariograma) é uma estatÃstica descritiva quantitativa que pode ser representada graficamente de tal forma que caracteriza a continuidade espacial de um conjunto de dados.
Em outras palavras, a análise variográfica caracteriza a (auto)correlação espacial.
Para esse cruso foi selecionado um estudo desenvolvido em áreas pertencentes à Fazenda de Ensino, Pesquisa e Extensão (FEPE), da Faculdade de Engenharia de Ilha Solteira (FEIS - UNESP), localizada no municÃpio de SelvÃria, estado do Mato Grosso do Sul.
No seu histórico, as áreas eram originalmente cobertas por vegetação nativa do Cerrado até a década de 1970, quando, em 1978, foram desmatadas e passaram a ser usadas para culturas anuais, como milho, soja, algodão e adubos verdes, até 1986. Durante os anos de 1986-1987, as áreas foram convertidas para os seguintes usos: floresta plantada de eucalipto (EU) e sistema silvipastoril (SI), uma floresta plantada de aroeira-vermelha (Myracrodruon urundeuva) em consórcio com capim braquiária (Brachiaria decumbens). A área de EU (Eucalyptus camaldulensis) foi formada em 26 de abril de 1986 e a área de SI foi formada em dezembro de 1987.
No momento da realização das avaliações de emissão de CO2 do solo (FCO2) as áreas já haviam passado por mais de \(30\) anos de mudança de conversão. As determinações foram realizadas no perÃodo de 03 de fevereiro a 17 de junho de 2017, Na imagem abaixo são apresentados os dias de avaliação e os valores de precipitação (chuva em mm) ocorridas no perÃodo.
Ao final do perÃodo de determinação da emissão de CO2 do solo, as amostras de solo foram coletadas e todos os atributos fÃsicos e quÃmicos foram determinados.
Aspectos gerais das áreas de estudo, silvipastoril (SI) e eucalipto (EU), ano de 2017, MunicÃpio de SelvÃria, Mato Grosso do Sul.
Um estudo clássico de geoestatÃstica normalmente começa com uma descrição detalhada da análise exploratória dos dados (mapas espaciais, histogramas). Para isso vamos realizas análise descritiva da variável emissão de CO2 do solo na área de eucalipto. (OLIVEIRA, 2018 )
# Instalar pacotes necessários
# install.packages("tidyverse")
# install.packages("sp")
# install.packages("gstat")
# install.packages("corrplot")
# install.packages("skimr")
# install.packages("agricolae")
# Carregar os pacotes
library(tidyverse)
library(sp)
library(gstat)
library(corrplot)
source("R/my-functions.R")
geo_fco2.rds disponÃvel na
pasta data. Utilize as funções glimpse e
skim, dos pacotes {dplyr} e
{skimr} para a apresentação de um resumo inicial do banco
de dados.# Lendo o banco de dados
dados_geo <- read_rds("data/geo_fco2.rds")
glimpse(dados_geo)
skimr::skim(dados_geo)
tratamento ?# Número de categorias de tratamentos
dados_geo |>
pull(tratamento) |>
unique()
# Filtrar as observações para o tratamento "EU"
dados_geo |>
filter(tratamento == "EU")
dados_geo |>
filter(tratamento == "EU") |>
ggplot(aes(x=x, y=y)) +
geom_point() +
theme_bw()
Observe a disposição dos pontos amostrais.
dados_geo |>
filter(tratamento == "EU") |>
mutate(
fco2_class = cut(fco2,4)
) |>
ggplot(aes(x=x, y=y, color=fco2_class, size=fco2_class)) +
geom_point() +
theme_bw() +
scale_color_viridis_d()
# Mostrar os dois gradeados....
dados_geo |>
mutate(
fco2_class = cut(fco2,4)
) |>
ggplot(aes(x=x, y=y,
color=fco2_class,
size=fco2_class)) +
geom_point() +
theme_bw() +
scale_color_viridis_d() +
facet_wrap(~tratamento)
N, Média, Mediana, MÃnimo, Máximo, Variância, Desvio Padrão, Assimetria, Curtose e Coeficiente de Variação.# EstatÃsticas Descritivas para fco2
dados_geo |>
filter(tratamento == "EU") |>
summarise(
N = n(),
Media = mean(fco2, na.rm = TRUE),
Mediana = median(fco2, na.rm = TRUE),
Menor = min(fco2, na.rm = TRUE),
Maior = max(fco2, na.rm = TRUE),
Variancia = var(fco2, na.rm = TRUE),
Desvio_Padrao = sd(fco2, na.rm = TRUE),
Assimetria = agricolae::skewness(fco2),
Curtose = agricolae::kurtosis(fco2),
CV = 100*Desvio_Padrao/Media
)
summarise e across.dados_geo |>
filter(tratamento == "EU") |>
summarise(
across(
fco2:k, my_estat_desc
)
)
A análise geoestatÃstica é baseada na teoria das variáveis regionalizadas, que é uma função numérica com distribuição espacial, que varia de um ponto a outro com continuidade aparente, mas cujas variações não podem ser representadas por uma função matemática simples (MATHERON, 1963).
Uma variável regionalizada é uma variável aleatória que assume diferentes valores, de acordo com a sua posição na área de estudo.
Se todos os valores de uma variável regionalizada forem considerados em todos os pontos dentro de uma área amostral, a variável regionalizada é apenas uma de infinitas variáveis aleatórias.
Esse conjunto é chamado de função aleatória (FA) e é simbolizado por \(Z(x_i)\). Na prática, quando retiramos uma amostra de solo em um local com coordenadas definidas, temos apenas uma única realização da função aleatória (FA).
Para estimar valores em locais não amostrados, deve-se introduzir as restrições de estacionaridade estatÃstica. A existência de estacionaridade permite que o experimento possa ser repetido mesmo que as amostras sejam coletadas em pontos diferentes, pois elas pertencem à mesma população, com os mesmos momentos estatÃsticos (VIEIRA, 2000).
Em resumo, temos que os métodos geoestatÃsticos são ótimos quando os dados são:
normalmente distribuÃdos e estacionários (média e variância não variam significativamente no espaço).
Desvios significativos da normalidade e da estacionariedade podem causar problemas, portanto é sempre importante começar a análise estudando o histograma ou algum gráfico similar para checar a normalidade e o mapa dos valores no espaço para checar a ocorrência de tendências significativas.
fco2. Adicione os
valores de média, mediana,
primeiro e terceiro quartil no histograma na
forma de linhas verticais.dados_geo |>
filter(tratamento == "EU") |>
pull(fco2) |>
quantile() # median() mean()
dados_geo |>
filter(tratamento == "EU") |>
ggplot(aes(x=fco2)) +
geom_histogram(bins=15,fill="lightgray",color="black") +
theme_bw() +
geom_vline(xintercept = c(4.65, 4.53, 3.67, 5.36),
linetype=2,
color= c("red","blue" ,"blue" ,"blue"),
lwd = 1)
dados_geo |>
filter(tratamento == "EU") |>
ggplot(aes(sample = fco2)) +
stat_qq() +
stat_qq_line(color="blue",lwd=2) +
theme_bw()
fco2.dados_geo |>
filter(tratamento == "EU") |>
pull(fco2) |>
shapiro.test() # teste de normalidade
x e y.dados_geo |>
filter(tratamento == "EU") |>
ggplot(aes(x=x, y=fco2)) +
geom_point() +
theme_bw()
dados_geo |>
filter(tratamento == "EU") |>
ggplot(aes(x=y, y=fco2)) +
geom_point() +
theme_bw()
dados_geo |>
filter(tratamento == "EU") |>
lm(formula = fco2 ~ x) |>
summary.lm()
dados_geo |>
filter(tratamento == "EU") |>
lm(formula = fco2 ~ y) |>
summary.lm()
As medidas de dispersão conjunta de duas variáveis, como X e Y, referem-se à extensão em que os valores dessas variáveis se afastam de seus valores médios simultaneamente. Algumas medidas comuns incluem a covariância e a correlação.
A covariância \((Cov)\) indica a direção da relação linear entre duas variáveis. Se for positiva, sugere uma relação positiva (quando uma variável aumenta, a outra também tende a aumentar); se for negativa, indica uma relação negativa (uma variável aumenta enquanto a outra diminui).
Onde \(X_i\) e \(Y_i\) são os valores individuais das variáveis \(X\) e \(Y\), \(\bar{X}\) e \(\bar{Y}\) são as médias de \(X\) e \(Y\), e \(n\) é o número de observações.
\[ Cov(X,Y) = \frac{1}{n}\sum\limits_{i=1}^n(X_i - \bar{X})(Y_i - \bar{Y}) \]
A correlação \((r)\) normaliza a covariância para estar entre \(-1\) e \(1\), facilitando a interpretação. Um valor próximo de \(1\) indica uma forte relação positiva, \(-1\) indica uma forte relação negativa, e \(0\) indica ausência de relação linear.
\[ r=\frac{\sum\limits_{i=1}^n(X_i - \bar{X})(Y_i - \bar{Y})}{\sqrt{\sum\limits_{i=1}^n(X_i - \bar{X})^2\sum\limits_{i=1}^n(Y_i - \bar{Y})^2}} \]
dados_geo |>
filter(tratamento == "EU") |>
select(fco2:k) |>
plot()
dados_geo |>
filter(tratamento == "EU") |>
select(fco2:k) |>
cor() |>
corrplot(method = "ellipse",
type = "upper")
dados_geo |>
filter(tratamento == "EU") |>
select(fco2:k) |>
GGally::ggpairs() +
theme_bw()
GeoestatÃstica lida com dados autocorrelacionados espacialmente.
Autocorrelação: correlação entre elementos de uma série com outros elementos da mesma série separados por uma uma determinada distância. (Oxford American Dictionary)
Covariância e a correlação são medidas de similaridade entre duas variáveis. Para estender este conceito para medidas de similaridade espacial, considere um gráfico de dispersão em que os pares de dados representam medidas da mesma variável separadas por uma distância \(h\) uma da outra.
Esta distância é geralmente denominada de lag, como
usado nas series temporais.
Ao contrário das funções de covariância e correlação, as quais medem a similaridade, o semivariograma experimental \(\gamma(h)\) mede dissimilaridade média entre dados separados por um vetor \(h\), ou seja,
\[ \gamma(h) = \frac{1}{2N(h)} \sum \limits_{i=1}^{N(h)}\left[ z(x_i) -z(x_i + h)\right]^2 \]
Sendo \(z(x_i)\) e \(z(x_i + h)\) os valores de cauda e da cabeça, respectivamente e \(N(h)\) o número de pares de pontos separados pela distância \(h\).
x, y e
z.x <- dados_geo |>
filter(tratamento == "EU") |>
pull(x)
y <- dados_geo |>
filter(tratamento == "EU") |>
pull(y)
z <- dados_geo |>
filter(tratamento == "EU") |>
pull(fco2)
matriz_dist <- matrix(0,ncol=102,nrow=102)
matriz_gamma <- matrix(0,ncol=102,nrow=102)
\[D = \sqrt{(X_2-X_1)^2+(Y_2-Y_1)^2}\] e
\[\gamma = \left[ z(x_i) -z(x_i + h)\right]^2\]
, respectivamente:
for(i in 1:102){
for(j in 1:102){
matriz_dist[i,j] = sqrt((x[i]-x[j])^2 + (y[i]-y[j])^2)
matriz_gamma[i,j] = (z[i]-z[j])^2
}
}
table(matriz_dist) %/% 2 |>
plot()
abline(h=30, col="red", lty =2)
recomenda-se um núemro de pares de ponto em uma dada distância de separação \(h\) sempre maior que \(30\).
lag <- 5
active_lag_distance <- 60
h <- seq(5, active_lag_distance, lag)
n_h <- h
gamma_h <- h
for(i in seq_along(h)){
if(i == 1){
filtro <- matriz_dist > 0 & matriz_dist <= h[i]
}else{
filtro <- matriz_dist > h[i-1] & matriz_dist <= h[i]
}
n_h[i] <- sum(filtro)/2
gamma_h [i] <- sum(matriz_gamma[filtro])/2/n_h[i]/2
}
tibble(n_h,h,gamma_h) |>
ggplot(aes(x=h,y=gamma_h)) +
geom_point(size=3) +
theme_bw()
O gráfico de acima é denominado semivariograma experimental, o qual exibe um comportamento puramente aleatório ou sistemático que pode ser descrito por modelos matemáticos teóricos (linear, esférico, exponencial, gaussiano e lei de potência).
O gráfico do semivariograma experimental, \(\hat{\gamma}(h)\), é formado por uma série de valores, sobre os quais se objetiva ajustar uma função. É importante que o modelo ajustado represente a tendência de \(\hat{\gamma}(h)\) em relação a \(h\).
O procedimento de ajuste não é direto e automático, como no caso de uma regressão, por exemplo, mas sim interativo, pois nesse processo faz-se um primeiro ajuste e verifica-se a adequação do modelo teórico. Dependendo do ajuste obtido, pode-se ou não redefinir o modelo, até obter um que seja considerado satisfatório.
Os principais modelos a serem ajustados são:
\[\begin{cases} \hat{\gamma}(h) = 0 \text{ se } h=0 \\ \hat{\gamma}(h) = Cte \text{ caso contrário } \end{cases}\]
$$ (h) = b(h)^c 0 < c < 2 h
$$
\[\hat{\gamma}(h) = C_0 + C_1 \left[1- exp\left({-3 \frac{h}{a}}\right) \right], h > 0\]
\[\begin{cases} \hat{\gamma}(h) = C_0 + C_1 \left[exp\left(\frac{3}{2}(\frac{h}{a})-\frac{1}{2}(\frac{h}{a})^3\right)\right], 0 \leq h \leq a \\ \hat{\gamma}(h) = C_0 + C_1, h>a \end{cases}\]
\[\hat{\gamma}(h) = C_0 + C_1 \left[1- exp(- \frac{h^2}{a^2}) \right], h > 0\] O modelo esférico alcança o valor do patamar \((C_0+C_1)\) no valor de alcance \((a)\) especificado. Os modelos exponencial e gaussiano se aproximam do patamar assintoticamente, com a representando o alcance prático, ou seja, a distância na qual a semivariância alcança \(95\%\) do valor do patamar.
Para variáveis que apresentam dependência espacial, espera-se que os valores de \(Z(x_i) – Z(x_i + h)\) aumentem com a distância \(h\) até determinada distância, a partir da qual os valores estabilizam.
O valor de semivariância, na qual ocorre a estabilização do semivariograma, é denominado de patamar (sill), representado pelo sÃmbolo \(C_0 + C_1\), sendo aproximadamente igual ao valor da estimativa da variância dos dados analisados.
A distância \(h\) na qual ocorre a estabilização do semivariograma é denominada de alcance (range), simbolizado por \(a\), e define o limite da dependência espacial (autocorrelação). Na prática, pontos cuja distância entre eles seja maior que o valor de \(a\), não apresentam autocorrelação espacial, ou seja, são aleatórios.
O efeito pepita (nugget effect), representado pelo sÃmbolo \(C_0\), é o valor de semivariância encontrado no intercepto do modelo ajustado com o eixo \(Y\). Teoricamente, este valor deve ser zero para uma distância de separação \((h)\) igual a zero; entretanto, erros de amostragem e a variabilidade na pequena escala podem causar desvio do zero para esse parâmetro. Portanto, o efeito pepita representa a quantidade de variância não explicada ou modelada como correlação espacial.
O valor \(C_1\) representa a estrutura de variabilidade espacial dos dados.
Modelo Exponencial - descreve uma \(FA\) que é bem errática à pequenas distâncias, e a sua correspondente realização apresentada na figura a seguir, evidencia considerável variabilidade na pequena escala.
Modelos Esférico - descreve uma \(FA\) que é menos errática em pequenas distâncias, e sua correspondente realização pode ser verificada abaixo.
Modelo Gaussiano descreve uma \(FA\) que é extremamente contÃnua, o seu semivariograma é tangencial ao eixo x na origem, e cresce suavemente. A sua correspondente realização, figura abaixo é suave e levemente ondulada.
coordinates() e
variogram() dos pacotes {sp} e
{gstat}, construa o semivariograma experimental \(h\) vs \(\gamma(h)\), não esqueça de adicionar o
número de pares de pontos de cada valor estimado de semivariância para
futuras inspeções.# Criando o arquivo necessário somente x, y, e z
df_aux <- dados_geo |>
filter(tratamento == "EU") |>
select(x, y, fco2) |>
rename(z=fco2)
class(df_aux)
coordinates(df_aux) = ~x + y
class(df_aux)
# Criando a fórmula para
formula <- z~1
# Criando o semivariograma experimental
semivariograma_experimental <- variogram(formula, data=df_aux,
width = 5,
cutoff=50,
cressie=TRUE)
semivariograma_experimental |>
ggplot(aes(x=dist,y=gamma)) +
geom_point() +
theme_bw() +
annotate("text", x=semivariograma_experimental$dist,
y=semivariograma_experimental$gamma*.98,
label = semivariograma_experimental$np)
fit.variogram, vgm e
variogramLine, pertencentes ao pacote
{gstat}.modelo_ajustado <- fit.variogram(semivariograma_experimental,
vgm(model = "Sph",nugget = 0.5,
psill = 1.2, range = 25 ))
preditos <- variogramLine(object = modelo_ajustado, maxdist = 50)
my_fitted_semivar <- semivariograma_experimental |>
ggplot(aes(x=dist,y=gamma)) +
geom_point() +
geom_line(data = preditos, aes(x = dist, y = gamma)) +
theme_bw()
my_fitted_semivar
# Extraindo os parâmetros do modelo
sse <- attr(modelo_ajustado, "SSErr")## Soma de Quadrado dos ResÃduos
c0 <- modelo_ajustado$psill[[1]]## Valor de Nugget
patamar <- sum(modelo_ajustado$psill)## Valor de Patamar
alcance <- modelo_ajustado$range[[2]]## Alcance
## Nomes dos parâmetros
nomes <- c("C0", "C0+C1", "Range", "SSE")
## vetor de valores
vetor_par <- c(c0, patamar, alcance, sse)
## Plotando o semivariograma
my_fitted_semivar +
annotate("text",x=rep(50,4),
y=seq(.8,1,.06),
label=paste0(nomes,"= ",round(vetor_par,4)),size=5) +
labs(x="h",
y = expression(paste(hat(gamma),"(h)")))
A ideia marcante da geoestatÃstica é bem simples. Ela consiste nos seguintes passos:
Defina uma área/local \(A\), considerada homogênea o suficiente para a garantir a interpolação dentro dela, ou seja,
Examine todos os dados medidos dentro de \(A\) para calcular as caracterÃsticas- \(h\) da variabilidade espacial, isto é, calcule os valores do semivariograma experimental.
Modele o semivariograma experimental com uma função avaliável para todos os vetores de distância \(h\).
Os passos mais importantes e, por isto mesmo, os que mais trazem conseqüências em qualquer estudo geoestatÃstico são:
A decisão de estacionariedade está implÃcita em todas as estatÃsticas, ela não é particular da abordagem geoestatÃstica. Esta decisão permite definir um conjunto de dados (área A) sobre os quais as médias experimentais e proporções serão calculadas e é assumido representativa da população como um todo, não de qualquer particular locação individual \((x \in A)\).
A hipótese intrÃnseca assume que o incremento
\[ Z(x+h) - Z(x) \]
tem média zero, ou seja,
\[ E[Z(x+h) - Z(x)] = 0 \]
e uma variância finita que não depende do local \(x\),
\[ Var[Z(x+h) - Z(x)] = E[Z(x+h) - Z(x)]^2 = 2\gamma(h) \]
As convenções direcionais usadas na geoestatÃstica são mostradas na Figura abaixo
MODELO ISOTRÓPICO - Semivariograma só depende da distância \(h\). Um único modelo modela a variabilidade em todas as direções.
MODELO ANISOTRÓPICO – A modelo aplicado ao semivariograma experimental depende da distância de separação e da direção.
O variograma ao longo da direção 30º, por exemplo, de um canal de sedimentação, apresentará menor variabilidade do que na direção 120º, perpendicular à continuidade, correspondentemente, é observado maior correlação na direção da sedimentação.
Estimativa em locais não amostrados
Use este modelo \(\hat{\gamma}(h)\), ou o modelo de covariância \(C(h)\), para a tradicional técnica de regressão, de krigagem.
Conhecido o semivariograma, e havendo dependência espacial entre as amostras pode-se estimar valores da propriedade \(Z\) em qualquer ponto do campo de estudo, sem tendência e com variância mÃnima. O método de estimação é denominado de KRIGAGEM, nome dado por MATHERON (1963) em homenagem ao engenheiro de minas sul-africano KRIGE.
### Estimação Local: Krigagem
\[ \hat{Z}(X_0) = \lambda_1Z(X_1) + \lambda_2Z(X_2)+ \cdots + \lambda_kZ(X_k) = \sum \limits_{i = 1}^k \lambda_iZ(X_i) \]
Os pesos são escolhidos de forma que a estimativa \(\hat{Z}(X_0)\) seja não viciada e a variância do erro \(\sigma_e^2\) da estimativa seja menor que qualquer outra combinação linear dos valores observados, ou seja:
\[ E[\hat{Z}(X_0) - Z(X_0)] = 0 \text{ e } \\ \sigma^2_e = Var[\hat{Z}(X_0) - Z(X_0)] = E[\hat{Z}(X_0) - Z(X_0)] \text{ seja mÃnima} \]
Trabalhando algebricamente estas duas expressões, chegamos às expressões que forma o sistema de krigagem, que pode ser esrito em função do semivariograma \(\gamma(h)\), sabendo-se que \(\gamma(h) = C(0) - C(h)\)
\[ \begin{cases} \sum \limits_i \lambda_i(x_i,x_j) - \mu = \gamma(x_i,x_0) i = 1,2, \cdots, N \\ \sum \limits_i \lambda_i = 1 \end{cases} \]
E a variância da estimação fica:
\[ \sigma^2_{\hat{Z}(X_0)} = \sum \limits_i^N \lambda_i\gamma(X_i,X_0) + \mu \]
Sistema de krigagem na notação matricial.
Matriz de semivariância entre as amostras \[[\gamma]=\begin{bmatrix} \gamma(X_1,X_1) & \gamma(X_1,X_2) & \cdots & \gamma(X_1,X_N) & 1 \\ \gamma(X_2,X_1) & \gamma(X_2,X_2) & &\gamma(X_2,X_N) & 1 \\ \vdots & & \ddots & \vdots & 1\\ \gamma(X_N,X_1) & \gamma(X_N,X_2) & \cdots & \gamma(X_N,X_N) & 1 \\ 1 & 1 & 1 & 1 &0 \\ \end{bmatrix}\]
Matriz de semivariância entre as amostras e o ponto a estimar
\[[b]=\begin{bmatrix} \gamma(X_1,X_0) \\ \gamma(X_2,X_0) \\ \vdots \\ \gamma(X_N,X_0) \\ 1 \end{bmatrix}\]
Matriz de ponderadores
\[[\lambda]=\begin{bmatrix} \lambda_1 \\ \lambda_2 \\ \vdots \\ \lambda_N \\ \mu \end{bmatrix}\]
a Relação entre os sistema é escrita como:
\[ [\gamma][\lambda]=[b] \text{ ou seja: } \\ \]
\[ \begin{bmatrix} \gamma(X_1,X_1) & \gamma(X_1,X_2) & \cdots & \gamma(X_1,X_N) & 1 \\ \gamma(X_2,X_1) & \gamma(X_2,X_2) & &\gamma(X_2,X_N) & 1 \\ \cdots & & \ddots & \cdots & 1\\ \gamma(X_N,X_1) & \gamma(X_N,X_2) & \cdots & \gamma(X_N,X_N) & 1 \\ 1 & 1 & 1 & 1 &0 \\ \end{bmatrix} \begin{bmatrix} \lambda_1 \\ \lambda_2 \\ \cdots \\ \lambda_N \\ \mu \end{bmatrix}= \begin{bmatrix} \gamma(X_1,X_0) \\ \gamma(X_2,X_0) \\ \cdots \\ \gamma(X_N,X_0) \\ 1 \end{bmatrix}\]
Então,
\[[\lambda] = [\gamma]^{-1}[b]\]
cuja solunção única é \[\hat{Z}(X_0) = \sum \limits_i^N \hat{\lambda}_iZ(X_i)\]
e a expressão da variância da forma matricial fica:
\[ \sigma_e^2=[\gamma]^t[b] \]
# Estrair os vetores das coordenadas
x <- dados_geo |>
filter(tratamento == "EU") |>
pull(x)
y <- dados_geo |>
filter(tratamento == "EU") |>
pull(y)
# definir a menor qual será a distância de separação entre os pontos
distancia <- 1 # 2 m
gradeado_adensado <- expand.grid(
x = seq(min(x),max(x),distancia),
y = seq(min(y),max(y),distancia)
)
# Plotar o gradeado amostral e o adensamento
dados_geo |>
filter(tratamento == "EU") |>
ggplot(aes(x=x, y=y, color="blue")) +
geom_point() +
geom_point(data=gradeado_adensado,
aes(x=x,y=y),color = "blue",
shape=3,size=.2) +
theme_bw()+
theme(legend.position = "none")
# Transforma o arquvo do gradeado adensado em um objt SpatialData
coordinates(gradeado_adensado) = ~ x + y
class(gradeado_adensado)
# remotes::install_github("https://github.com/r-spatial/gstat")
ko <- krige(formula, df_aux, gradeado_adensado, modelo_ajustado,
block =c(0,0),
nsim = 0,
na.action = na.pass,
debug.level = 1)
as_tibble(ko) |>
ggplot(aes(x,y)) +
geom_tile(aes(fill=var1.pred)) +
scale_fill_viridis_c()
as_tibble(ko) |>
ggplot(aes(x,y)) +
geom_tile(aes(fill=var1.var)) +
scale_fill_gradient(low = "yellow",high = "blue")
GOTWAY, C.A.; HARTFORD, A.H. 1996. Geostatistical methods for incorporating auxiliary information in the prediction of spatial variables. J. Agric., Biol. Environ. Statis., 1: 17-39.
GOOVAERTS, P. Geostatistics for natural resources evaluation. New York: Oxford University Press, p. 483, 1997.
ISAAKS, E. H.; SRIVASTAVA, R. M. Applied geostatistics. Nova York: Oxford University Press, 1989. 561 p.
JOURNEL, A.G. Fundamentals of geostatistics in five lessons. Short Course in Geology: Volume 8. American Geophysical Union, Washington, p. 1 – 40, 1989.
MATHERON G. Principles of geostatistics. Economic Geology, 58, 1963, 1246-1266.
OLIVEIRA, C. F. Variabilidade espacial da emissão de CO2 e estoque de carbono do solo em áreas de eucalipto e sistema silvipastoril. Dissertação de Mestrado, Agronomia, UNESP-FEIS, 2018.
VIEIRA, S. R. GeoestatÃstica em estudos de variabilidade espacial do solo. In: NOVAIS, R.F.; ALVAREZ, V. H.; SEHAFFER, C. E. G. R. (Ed). Tópicos em ciência do solo. Viçosa: Sociedade Brasileira de Ciência do Solo, 2000, v.1, p. 1 – 53.
WEBSTER, R., OLIVER, M. Geostatistics for Environmental Sciences. San Francisco: John Wiley & Sons, p. 315, 2009.
A todos os participantes e integrantes do grupo de pesquisa.